##################################################
#
# PLOTS OF HARSH REPLICATION MODELS FOR
# OIL, ISLAM, and WOMEN PROJECT
#
# Paul Musgrave and Yu-Ming Liou
# musgrave@umass.edu and yl254@georgetown.edu
#
# Created 27 February 2016
#
##################################################


##################################################
# R Housekeeping
##################################################

# Remove and empty workspace
rm(list=ls())

#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 
# Load Packages
#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 

library(foreign)
library(sandwich)
library(apsrtable)
library(tonymisc)  # includes summaryR()
library(MASS) # for polr() which runs our ologit

#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 
# Load Functions
#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 


# Change this to  your working directory
mywd <- "/Users/paulmusgrave/Dropbox/0001 Academic Projects/Ongoing/0127 Oil Islam Women/ISQ Accepted Submission/Replication"

setwd(paste(mywd,"/Code",sep=""))

source("statifyFN.R")
source("multipleOLSFN.R")
source("genModelsFN.R")
source("plotlineFN.R")
source("plotprepFN.R")
source("modelDefinitions.R")

#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 
# Load Data
#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 

setwd(paste(mywd,"/Data",sep=""))

data <- read.dta("NewDVModelsXSforISQFINAL2015.dta")
data$logGDPcap6372_sq <- data$logGDPcap6372 * data$logGDPcap6372

#################################################
#
# Models of Oil Income and Reproductive Issues
# 
#################################################

# Use p140 alt

out.totfert <- multipleOLS(formulas=genModels(c("totfert"),
                                            c(p140alt,p142base)),
				          	modelnames=c(names(p140base),names(p142base)))
# Note modelnames changed for this and others
out.adfert <- multipleOLS(formulas=genModels(c("adfert"),
                                             c(p140alt,p142base)),
                          modelnames=c(names(p140base),names(p142base)))

out.contraception <- multipleOLS(formulas=genModels(c("avgcontraception"),
                                             c(p140alt,p142base)),
                          modelnames=c(names(p140base),names(p142base)))


## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## #
# Plotting Models
## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## #
vars.to.plot <- c("oil_gas_valuePOPred")
				
plot.out.totfert <- do.call(rbind.data.frame,
                            lapply(out.totfert,plotprep,vars.to.plot))
plot.out.adfert <- do.call(rbind.data.frame,
                           lapply(out.adfert,plotprep,vars.to.plot))
plot.out.contraception <- do.call(rbind.data.frame,
                            lapply(out.contraception,plotprep,vars.to.plot))

plot.out.totfert$Index <- 1:nrow(plot.out.totfert)
plot.out.adfert$Index <- 1:nrow(plot.out.adfert)
plot.out.contraception$Index <- 1:nrow(plot.out.contraception)

plot.out.contraception$bg <- ifelse(abs(plot.out.contraception$Estimate/plot.out.contraception$SE)>=1.96,"black",
                                        ifelse(abs(plot.out.contraception$Estimate/plot.out.contraception$SE)<=1.65,"white",
                                               "gray"))    

plot.out.adfert$bg <- ifelse(abs(plot.out.adfert$Estimate/plot.out.adfert$SE)>=1.96,"black",
                                    ifelse(abs(plot.out.adfert$Estimate/plot.out.adfert$SE)<=1.65,"white",
                                           "gray"))    
plot.out.totfert$bg <- ifelse(abs(plot.out.totfert$Estimate/plot.out.totfert$SE)>=1.96,"black",
                             ifelse(abs(plot.out.totfert$Estimate/plot.out.totfert$SE)<=1.65,"white",
                                    "gray"))    
# pch = 21
# col = black

setwd("../Drafts/Charts")

# FOR PDF
#pdf(file="ISQFINALMainBodyFigure4",height=5,width=7)

# FOR HI-RES TIFF
tiff(file="ISQFINALMainBodyFigure4.tiff",
		height=5,
		width=7,
		units="in",
		res=300)

par(mfrow=c(3,1),oma=c(0,0,3,1),mar=c(3,3,1,1))

plot(1,
     pch="",
     ylim=c(-1.15*min(plot.out.totfert$Y0),1.15*max(plot.out.totfert$Y1)), 
     xlim=c(.5,max(plot.out.totfert$Index)+.5), 
     xaxt="n",
     xlab="",
     ylab="Coefficient Estimate",
     main="Estimated Relationship With Total Fertility",
     cex.main=1.25,
     adj=1)
abline(h=0,lty="dotdash")

segments(y0=plot.out.totfert$Y0,
		y1=plot.out.totfert$Y1,
		x0=plot.out.totfert$Index,
		x1=plot.out.totfert$Index,
		col="black",
		lwd=1,
		lty=1)
points(x=plot.out.totfert$Index,
		y=plot.out.totfert$Estimate,
		col="black",
		pch=21,
		bg=plot.out.totfert$bg,
		cex=1.5)
axis(side=1,at=plot.out.totfert$Index,
		labels=plot.out.totfert[1:nrow(plot.out.totfert),]$Model,
		cex.axis=1.25, tick=TRUE)
		

### And Adolescent Fertility
plot(2,
     pch="",
     ylim=c(1.15*min(plot.out.adfert$Y0),1.15*max(plot.out.adfert$Y1)), 
     xlim=c(.5,max(plot.out.adfert$Index)+.5), 
     xaxt="n",
     xlab="",
     ylab="Coefficient Estimate",
     main="Estimated Relationship With Adolescent Fertility",
     cex.main=1.25,
     adj=1)
abline(h=0,lty="dotdash")

segments(y0=plot.out.adfert$Y0,
         y1=plot.out.adfert$Y1,
         x0=plot.out.adfert$Index,
         x1=plot.out.adfert$Index,
         col="black",
         lwd=1,
         lty=1)
points(x=plot.out.adfert$Index,
       y=plot.out.adfert$Estimate,
       col="black",
       pch=21,
       bg=plot.out.adfert$bg,
       cex=1.5)
axis(side=1,at=plot.out.adfert$Index,
     labels=plot.out.adfert[1:nrow(plot.out.adfert),]$Model,
     cex.axis=1.25, tick=TRUE)


plot(3,
     pch="",
     ylim=c(1.15*min(plot.out.contraception$Y0),
                abs(1.15*max(plot.out.contraception$Y1))), 
     xlim=c(.5,max(plot.out.contraception$Index)+.5), 
     xaxt="n",
     xlab="",
     ylab="Coefficient Estimate",
     main="Estimated Relationship With Contraceptive Availability",
     cex.main=1.25,
     adj=1)
abline(h=0,lty="dotdash")

segments(y0=plot.out.contraception$Y0,
         y1=plot.out.contraception$Y1,
         x0=plot.out.contraception$Index,
         x1=plot.out.contraception$Index,
         col="black",
         lwd=1,
         lty=1)
points(x=plot.out.contraception$Index,
       y=plot.out.contraception$Estimate,
       col="black",
       pch=21,
       bg=plot.out.contraception$bg,
       cex=1.5)
axis(side=1,at=plot.out.contraception$Index,
     labels=plot.out.contraception[1:nrow(plot.out.contraception),]$Model,
     cex.axis=1.25, tick=TRUE)

# Title for whole thing
mtext("Estimating Relationship Between Oil Income And Reproductive Behavior", side=3, line=1, outer=TRUE, cex=1.15, font=2)

dev.off()
		
		
###